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COMPLETE GALILEAN-INVARIANT LATTICE BGK MODELS FOR THE 

NAVIER-STOKES EQUATION 

YUE-HONG QIAN* AND YE ZHOUt 

Abstract. Galilean invariance has been an important issue in lattice- based hydrodynamics models. 
Previous models concentrated on the nonlinear advection term. In this paper, we take into account the 
nonlinear response effect in a systematic way. Using the Chapman-Enskog expansion up to second order, 
complete Galilean invariant lattice BGK models in one dimension ( 0 = 3) and two dimensions (0 = 1) for 
the Navier-Stokes equation have been obtained. 

Key words. Boltzmann equation, lattice-based hydrodynamics models, Navier-Stokes equation, Galilean 
invariance 

Subject classification. Fluid Mechanics 

1. Introduction. The modern version of Lattice Gas Automata (LG A) hydrodynamics model, intro- 
duced by Frisch, Hasslacher and Pomeau [1-2], has attracted much attention since 1986. In particular, the 
natural advantage of lattice-based models for performing parallel computation and handling complex ge- 
ometries problems, with local operations, indicates the lattice models may provide an alternative efficient 
numerical scheme for studying many complex phenomena [3-8]. 

However, the disadvantages of LGA models are rather significant. Boolean operations are intrinsically 
noisy. The pressure also contains an explicit velocity dependence, which is obviously unphysical. Further- 
more, LGA could not satisfy Galilean invariant due to the density dependence of a so-called g{p) factor [9], 
We stress that Galilean invariance is a basic requirement that is demanded for any physical model. 

Steps were taken to overcome the difficulties of LGA. McNamara and Zanetti [10] replaced the Boolean 
operations by neglecting particle correlations and introducing averaged distribution functions. Therefore, 
the computational accuracy is improved by eliminating the statistical noise. Higuera and Jimenez [11] noted 
that there are many choices for the collision integral that could generate the same equilibrium distribution. 
In 1992, the Bhatnagar-Gross-Krook approximation for lattice Boltzmann was independently introduced by 
Qian et al. [12-13] and Chen et al. [14-15]. The key observations in these two studies are that the requirement 
for using Fermi-Dirac statistics is no longer necessary and that of the resulted freedom in the functional form 
of the equilibrium distribution provide many choices for using desired equilibrium distribution functions. 
The Galilean invariance can be easily achieved for advection terms. 

Qian and Orszag [16] found a cubic nonlinear deviation from the Navier-Stoke equations when nonlinear 
response was taken into account. The cubic nonlinear term is important in compressible regime when the 
Mach number is not small, as usually assumed in the LGA models. Note that the cubic term, as a direct 
result of popular equilibrium functions, is hidden in all lattice Boltzmann calculations. Note that Chen et 
al. [17] attempted to eliminate the cubic nonlinear deviation term. Their effort was motivated by the fact 
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that the cubic terms is not part of original Navier- Stokes equation. They used an artificial procedure by 
introducing and manipulating higher order tensors of the expanded equilibrium distribution function in a 
rather complicated manner. 

The purpose of this letter is to point out that current lattice Boltzmann and lattice BGK models have not 
achieved a complete Galilean invariance. Indeed, we must emphase here that the cubic nonlinear deviation 
term is not Galilean invariant when u — ► u-j-U 0 where U 0 is a constant velocity. We will present an expanded 
discussion on the Galilean invariance issue in next section and develop a 5- velocity model in one dimension 
( 0 = 3.0) and a 17- velocity model in two dimensions (9 = 1.0) to eliminate the cubic nonlinear term. As a 
result, these two models achieve a status of complete Galilean invariance up to second order derivative since 
no higher- order dynamics exists [18]. Along the way, these two models remove another shortcoming of the 
lattice Boltzmann schemes, the staggered spurious invariants, as pointed out by the paper [19]. 

2. Galilean Invariance and Models Description. The original FHP model [1] has several weak- 
nesses. Non-Galilean invariance is one of such problems. Much effort has been made to get the Galilean 
invariance in the nonlinear advection term, first by lattice gas models [9], then by lattice BGK models [12- 
15]. Qian [12] pointed out that the problem of non-Galilean invariance is directly related to the problem 
of velocity-dependent pressure, the two problems became one. When deriving large-scale dynamics, the 
models use linear response of small perturbation. If nonlinear response effect is taken into account, a term 
deviating from the Navier-Stokes equation is obtained [16]. This new term is non-Galilean invariant and 
this fact hasn’t been pointed out so far as we know. The consequence of this cubic nonlinear term leads to 
a frame- velocity dependent viscosity, which is not physical. We shall recall first the previous models before 
introducing a tunable parameter 6 to control the cubic nonlinear term. 

Time evolution of lattice BGK models is the following [13], 

( 2 . 1 ) fi{x + Ci,t+ 1 ) = fi(x,t) +w[f?(x,t) - fi{x,t)\ 


where /* is the particle density with given velocity cj and lj the relaxation parameter (0 < u < 2). i is the 
index of discrete velocity c*. A workable equilibrium ff leading to hydrodynamic equations is [12-13], 


( 2 . 2 ) 
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where t p is a lattice structure weighting factor (the index p = cf) and c a is a constant. A slightly different 
form of equilibrium distribution for the FHP hexagonal lattice was used by Chen et al. [15]. The conditions to 
determine the values of t p are the constraints of conservation laws, the isotropic fourth order tensor of velocity 
and Galilean invariance in advection term. Greek subscripts a,/? denote the space directions in Cartesian 
coordinates. The conventional summation of repeated Greek subscripts is used. The hydrodynamic density 
p and velocity u are defined as, 

( 2 . 3 ) p =Y^fi=Y,ft, pu = J 2 f& = Ylfi 5i 

i i i i 


Using the Chapman-Enskog expansion, we derive the hydrodynamic equations up to the second order of 
Knudsen number at long wavelength and long time limits, 


(2.4) d t p 4- d a (pu a ) = 0 

dt(pUa) + d 0 (pu a u/ 3 ) = -d a P + dp[pud 0 {u a ) + ptd a {up)} 

(2.5) —crdpd~ l (pu a U 0 U~ l ) 
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where v and £ are shear and bulk viscosities [13], 


(2.6) "“‘“f (!-‘) 

The equation of state P{p) is, 

(2.7) P = c\p 


u is directly related to cu by, 


( 2 . 8 ) 


a = 


l 

CJ 


1 

2 


The (J cubic nonlinear term is the new term in eq. (2.5) which has been largely ignored in the community 
[16], the physical meaning of this term is the nonlinear response due to the quadratic term in the equilibrium 
( 2 . 2 ). 

We argue that this term is also a source of non- Galilean invariance. With a uniformly moving frame 
velocity ?7o, we have the following variable transformation: t — > t, x — * x — Uot , u — > u ± Uq and p — ► p. The 
momentum in the moving frame is now, 


dt(pua) + dp(pu a up) = -d a P + dp\pvdp(u a ) + p(;d a (up)] 

(2.9) —adpdj[p{u + U 0 ) a {u + U 0 )p(u + U 0 )^\ 

When measuring the shear viscosity, we linearize the above equations by assuming that Uq is in x direction, 

(2.10) d t v = {v- aUl)d x d x v 


where v is the transverse small velocity. Therefore, the effective shear viscosity i/ e ff is, 

( 2 . 11 ) = 

The velocity-dependence of viscosity is also unphysical. We will verify this relation numerically. 

We now use a modified equilibrium distribution for a one-dimensional 5- velocity model (D1Q5): (0, ±1, ±2) 
and a two-dimensional 17-velocity model (D2Q17): (0,0), (±1,0), (0, ±1), (±1,±1), (±2,0), (0, ±2) and 
(±2, ±2). We add a cubic nonlinear term to the equilibrium distribution with an ajustable parameter 6 , 


( 2 . 12 ) 
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In order to determine the weighting factors t p , in addition to the previous constraints the isotropic 6 th order 
tensor of particle velocity (for two- and three-dimensional cases) must be satisfied. Therefore, the problem 
to eliminate the cubic nonlinear term in (2.5) becomes how to determine the weighting factors t p and the 
parameter 6 in addition to the usual constraints [12-13], 

The final momentum equation for D1Q5 model is, 


(2.13) d t {pu) + d x (pu 2 ) = - d x P + d x [p{u + <)cku] + X^x9 x (pu 3 ) 


where x 

(2-14) x={e Lh^4 _ l)a 
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Fig. 2.1. Effective Shear viscosity v in function of the frame-velocity Uo. The DSQ17 model with 0 = 0.0, 2.0 (non- 
Galilean invariant ) and with 0 = 1.0 (Galilean invariant) and uj = 1.0 are used. Squares are numerical results and solid curves 
are predictions. 


6 = 3.0 leads to \ = 0.0. 

For the D2Q17 model, 

dt(pUa) + dp{pu a up) = -d a P + dp[pvdp(u a ) + pCd a (u 0 )\ 

(2.15) +xdpd- l (pu a u 0 u 1 ) 

and x is also derived, 

(2.16) X = - 1)^ 

We need to replace — a in Eq. (2.10) with x and 6 = 1.0 leads to x = 0-0. After some algebraic manipulations, 
the c* for the D2Q17 model is the solution of the following quadratic equation, 

(2.17) hlX 2 + 13X — 12 = 0 

We obtain in the following table the weighting factors of these two models, 


Table 1 

Weighting factors t Pf constant c 8 and 0 


Model 

to 

t\ 

1 2 

<4 

* 8 


0 

D1Q5 

1/2 

1/6 


1/12 


i 

3.0 

D2Q17 

0.4092905 

0.1123018 

0.0335591 

0.0017273 

0.0000891 

0.3740845 

1.0 


Using these parameters, the cubic nonlinear term is exactly cancelled out without introducing any higher 
order nonlinearities. Therefore, the final equations are, 

(2.18) d t p + d a (pu a ) = 0 

(2.19) p(d t u a + updpu a ) = -d a P + dp\pvdp(u a ) + pC>da{up)\ 
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x-distance (dIQ5 Model with w-1, 0-0 and 0-3) 



Fig. 3.1. Comparison of shock tube simulations with (6 = 0.0^ or without (0 ~ 3.0) the cubic nonlinear term in the 
momentum equation. The D1Q5 model is used, lj = 1.0. 


which lead to the complete Galilean invariance regardless of the velocity direction of a uniformly moving 
frame. 

3. Simulations. Two numerical simulations have been carried out. The first one is to measure the 
effective shear viscosity in different moving frames. Figure 1 shows the effective shear viscosity in function 
of the uniform velocity Uq with 0 = 0.0, 1.0 and 2.0. Squares are numerical simulations while solid curves 
are predictions based on Eq. (2.10). Only 0 = 1.0 leads to velocity- independent viscosity which achieves the 
complete Galilean invariance. The neglect ion of this term results in relative errors of about 10% for Uq = 0.2. 
A shock tube simulation [20] is given in Figure 2 with the following conditions: initially, left half-tube density 
is 3.0 while right-half density 1.0, velocity is zero everywhere; at the two ends, bouncing-back conditions are 
imposed. The D1Q5 model is used to compare the cubic nonlinear effect with 0 = 0.0 and 0 = 3.0. The 
density and velocity differences are obtained at two different time steps t ~ 100 and t = 200. About 10% of 
difference is found due to the non- Galilean invariance. We do believe that non- Galilean invariant effect may 
be more significant in convection of passive scalars as discussed in the paper [9]. 

4. Conclusion. In this letter, we note that the current popular lattice-based hydrodynamics models 
are only Galilean invariant in the advection terms. The cubic nonlinear deviation term, which is a direct 
result of the chosen equilibrium distribution, does not satisfy the Galilean invariance. This fundamental 
property in required for all physical system. We develop a one-dimensional 5- velocity model ( 6 = 3.0) and 
a two-dimensional 17- velocity model (6 = 1.0) that restore the Galilean invariance by eliminating the cubic 
term. Therefore, the original Navier- Stokes equation is recovered with complete Galilean invariance. A 
three-dimensional model with complete Galilean invariance is under investigation and a lengthy paper of the 
details and applications of the present models will be published in the future [21]. An important by-product 
of this model is the elimination of the staggered spurious invariants [19]. For cases in which the continuum 
limits (long wavelength and long time) break down, effects such as rarefied gas effect [22] and finite lattice 
size effect [23] occur. 
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